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I . INTRODUCTION 


Many  aspects  of  terminal  ballistics  testing  are  in  apparent  need  of 

procedural  standardization.  We  attempt  here  to  deal  in  part  with  one 

such  area:  representation  of  the  relationship  between  striking  and 

residual  velocities  implicit  in  sets  of  (v  , v ) data. 

s r 

A form  is  proposed  as  being  sufficiently  simple  and  versatile  to 
usefully  and  realistically  model  velocity  interdependence.  A regression 
technique  (a  direct  non-linear  least  squares  algorithm)  is  presented 
which,  through  empirical  parameter  determination,  establishes  a systematic 
method  for  obtaining  an  explicit  functional  relationship;  in  particular, 
a "limit  velocity1'  is  thereby  routinely  generated  from  input  data.  It 
is  suggested  that  what  is  here  described  might  substantially  contribute 
towards  a standardized  methodology  for  assessment  and  categorization  of 
velocity  data  and  towards  standardization  of  strategy  in  testing- 

In  short,  if  the  overall  procedure  is  accepted,  then  a formula  is  at 
hand  to  concisely  contain  experimental  velocity  information.  Indeed 
much  of  the  information  from  a data  set  is  condensed  into  a triple  (a, 

»£,  p)  of  numbers  which  uniquely  specifies  the  relationship.  Ideally, 
and  in  the  nature  of  future  intent,  these  numbers  should  be  qualified 
according  to  data  variance  and  endowed  with  a statistical  appraisal  of 
confidence. 

The  practicable  essence  of  this  report  is  cast  in  a comprehensive 
list  of  Fortran  directives  which  methodically  induces  machines  to  extract 
parameters  from  data  and  to  plot  corresponding  (v  , v ) curves. 


II.  A PROPOSED  FORM 
1.  Context  and  Form  Considerations 


A frequent  test,  of  multiple  special  interest,  in  terminal  ballistics 
consists  of  the  firing  of  a number  of  nominally  identical  projectiles 
(penetrators)  into  as  many  nominally  identical  targets  with  all  controlled 
phenomena  except  for  striking  velocity*,  vs,  being  nominally  invariant. 

The  measured  response  to  each  such  impact  which  is  of  present  concern  is 
projectile  residual  velocity,  vr<  A collection  of  data  points  (vs,  v_) 
is  thereby  obtained.  In  accord  with  common  usage,  each  shot  is  deemed 
to  result  in  "penetration,"  while  "perforation"  is  signalled  by  the 
criterion  vr  > 0. 


For  present  purposes  we  regard  striking  velocity  as  a controlled  variable . 
Also,  in  deference  to  prevailing  custom,  velocities  are  taken  to  be 
non~negative  numbers  rather  them  vectors. 


Pivrodintr  mira  htanir 


S 


Let  us  suppose  a given  projectile-target  situation,  by  which  it  is 
assumed  that  a given  projectile  impacts  a given  target  with  variable 
striking  velocity  v but  with  all  other  relevant  pre-impact  characteris- 
tics being  given  ana  constant.  From  this  point  and  this  precise  the 
discussion  admits  of  mathematical  idealization.  In  the  given  context 
residual  velocity  can  be  considered  a function  of  striking  velocity; 
vr  = f (vs)  for  vs  > 0.  We  assume  the  existence  of  a limit  velocity,  v^, 
characterized  (i.e.,  defined)  by  the  properties:  v = 0 for  0 < v < v 

1*  S jC 

and  v » 0 for  v > vr  ; equivalently,  v.  = max  {v  :v  = 0}  = inf  {v  :v 
r sfn  7 * i.  sr  sr 

> 0).* 


As  further  characteristics  of  a suitable  (and  workable)  model,  we 
desire  continuity  and  that,  for  vs  > v^,  the  function  be  strictly  in- 
creasing, smooth  (differentiable)  and  concave**  (every  chord  lies  beneath 
its  intercepted  arc).  The  latter  propert.es  are  secured  by  requiring 
that  the  function  possess  at  each  point  on  its  support  ({x  : f(x)  T 0}) 
a positive  first  and  a negative  second  derivative.  Such  stipulation  is 
of  course  largely  contrivance  and  is  not  experimentally  decidable  but 
is  convenient,  suggestive  of  physical  experience,  and  imposing  of  no 
apparent  conflict  with  experiment  or  theory. 

2 . Further  Form  Considerat ions : Simple  Penetration  Theory  and  the 

Hyperbola 

Traditional  penetration  theory,  especially  where  an  explicit  model 
is  developed,  has  been  most  extensive  (and  evasive  of  much  unpleasantness) 
in  dealing  with  the  case  of  a rigid  (mass-preserving),  essentially  non- 
deforming penetrator;  e.g.,  a steel  projectile  impacting  a thin  aluminum 
plate.  Such  treatments,  otherwise  tolerating  a diversity  of  assumptions 
about  the  physics  of  the  particular  penetration  process  conjectured, 
have  been  numerous. 


Ic  is  a pertinent  observation,  and  is  perhaps  not  generally  recognized, 
that  the  penetration  models  evolved  from  these  theories  have  almost 
invariably  adhered  to  one  basic  form;  specifically: 


0,  0 < v 

< V 
s 

, 2 

2. 

a (vs  - 

r~\ 

1 

> 

■ si  f-'.  call 

y ex  a 

l 

1/2 


(1) 


a the  possibility  of  modeling 
a "shatter  regime,"  i.e.,  wc  ape  insisting  that  if  f(xp  > 0 and 
x9>  Xp  then  f(x2)  > 0.  Tnis  is  not  considered  a serious  compromise 
wzth  the  generality  for  which  we  strive. 


The  inf(imun)  of  a set  is  its  greatest  lower  bound, 

** 

This  precludes  the  possibility  of  inflection  points. 
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For  instance,  the  Ponce let-Morin  hypothesis*  of  1S29  concerning  the 
resistance  offered  a projectile  in  motion  through  a dense  medium  leads, 
with  supplementary  assumptions,  rather  easily  to  this  form  (c.f.  Robertson 
[1943]);  so  also  does  the  quite  different  simple  energy -momentum  analysis 
of  Recht  and  Ipson  [1963],  A variety  of  other  approaches  falls  into  line 
(Nishiwaki,  Thomson,  Zaid  and  Paul,  etc.)*  The  various  models  ultimately 
differ  precisely  in  so  far  as  do  the  formulations  for  a and  v^;  which 
quantities  we  will  tend  to  view  as  parameters  to  be  subjected  to  optimal 
adjustment  in  a given  situation. 


Our  examination  of  available  data  tends  to  confirm  that  experimental 
results  can  often  be  well  represented  within  the  framework  of  form  (1), 
particularly  in  situations  where  there  is  not  (or  is  not  expected  to  be) 
excessive  projectile  deformation.  Theory  and  experience  persuade  us 
that  in  many  instances  (1)  does  offer  a viable  model.  We  will  feel 
obliged  presently  to  adopt  many  characteristics  of  this  form,  and  the 
form  itself  as  a special  case,  in  a proposed  more  versatile  basic  form. 


We  note  that  form  (1)  meets  the  requisites  of  the  previous  section. 
Additionally,  there  is  an  infinite  right  derivate  at  v^;  i .e.  , while 


d v 


does  not  exist  at  v^,  it  is  true  that  vl^m 


d v 


s *7  d v 


Also  av_ 


approximates  for  large  values  of  vg ; precisely,  vl^.®„  (a  vs  - vr)  = 0* 


Observe  that  for  vg  i-  v^,  (1)  is  a quadrant  of  a hyperbola  with 
center  at  the  origin,  major  axis  coincident  with  the  horizontal  axis, 
having  as  one  asymptote  the  line  through  the  origin  of  slope  a. 

The  assumption  of  projectile  rigidity  typically  associated  with  this 
form  facilitates  dealing  in  terms  of  kinetic  energy;  within  the  perfora- 
tion regime,  (1)  implies  linearity  between  striking  and  residual  energies: 

e = of  (ec  - e/,  where  symbol  meanings  re  presumably  evident.  It  may 

at  this  point  be  worth  noting  that  even  in  cases  where  significant  de- 
formation and  mass  loss  are  experienced,  we  have  observed  tenative 
experimental  indication  (from  situations  where  residual  mass  measurements 
are  available)  of  approximate  linear  correlation  between  striking  and 
residual  energies. 


A nodtLficatior*  of  earlier  "sectional-pressure”  theories  asserting  that 
the  force  resisting  penetration  at  a given  target  depth  is  jointly 
proportional  to  a depth- dependent  cross-sectional  projectile  area  and 
to  a linear  function  of  instantaneous  projectile  energy . 
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Since  vy  cannot  exceed  vs,  an  imposed  parameter  constraint  will  be 
0 < a < 1.  Kn  the  event  of  perforation  by  a rigid  projectile,  an  impli- 
cation of  the  special  case  a * 1 is  that  energy  lost  to  the  penetration 
process  is  independent  of  striking  velocity.  The  various  proposed  models 
of  type  (1)  generally  concur  that  a should  be  near  one  for  thin  targets 
and,  other  things  being  equal,  should  decrease  with  increasing  target 
thickness.  For  example,  from  the  formulation  of  Recht  and  Ipson  for  an 
assumed  plugging  mode  of  failure,  we  have  a = (1  ♦ r)'1  where  r is  the 
ratio  of  the  mass  of  the  ejected  plug  (small  for  thin  plates)  to  that  of 
the  projectile. 

Otherwise  extant  simple  penetration  models  are  prone  towards  abject 
empiricism  - randomly  selected  theories  alleged  to  fit  experimental  data, 
limited  of  scope,  curiously  behaved,  and  in  any  event  not  sufficiently 
adaptable  for  our  purposes  to  be  usefully  exploited*. 

3.  Present  Proposal : The  Basic  Form 

There  is  clear  experimental  indication  that  in  general  a somewhat 
more  flexible  model  is  required  to  represent  the  velocity  relationship 
being  considered.  In  particular  it  has  been  found  that  form  (1)  is 
often  not  capable  of  adequately  reflecting  observed  behavior  for  vs  above 
and  in  the  vicinity  of  the  derived  v^;  such  deficiency  can  be  rather 
dramatically  apparent  for  instance  in  multiple  plate  situations  and  in  the 
case  of  long  rods  perforating  comparatively  thick  armor.  Nor  is  this 
especially  surprising  for  it  is  just  in  this  region  that  one  might 
suppose  deformation  to  be  most  consequential;  it  is  expected,  c.f. 
Defoumeaux  [1973],  that  the  proportional  contribution  of  plastic  defor- 
mation (and  of  friction)  to  total  energy  expenditure  diaijushes  towards 
zero  as  vs  becomes  large.  In  addition  to  the  vagaries  introduced  by  the 
possibility  of  projectile  deformation,  the  general  inadequacy  of  (1) 
might  be  reflected  as  well  in  variability  of  target  deformation.  More 
significant  than  the  absolute  susceptibility  to  deformation  of  either 
projectile  or  target,  we  expect,  may  be  the  manner  in  which  the  defor- 
mation patterns  are  distributed  with  striking  velocity. 

A basic  form  assimilating  desired  general  characteristics**  and 
considered,  from  a wide  spectrum  of  experimental  evidence,  well  suited 
to  effectively  represent  observed  behavior  is: 


Some  frequently  encountered  phenomna  rf  the  no  models  which  we  find 
restrictive  or  at  variance  with  our  physical  expectations  are  unaccount- 
able discontinuities  and  inflection  points , either  lack  of  an  asymptote 
or  insistence  that  the  asymptote  have  always  a slope  of  one  (or  some 
other  immutable  value  chosen  by  caprice ) end  absence  of  an  infinite 
right  slope  (derivative)  at  \f. 

A partial  exception  will  be  noted  in  the  next  section. 
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(2) 


0,  0 < vs  < 


i/r 

* Cv/  - v£-)  . vs  > vz 


p _ p 

with  constraints  0 < <2  < 1 and  p > 1*. 


It  is  regarded  that  (2)  offers  additional  and  sufficient  versatility 
precisely  in  the  region  where  (1),  which  is  a special  case  of  (2),  was 
found  to  be  deficient.  The  parameter  a assumes  its  previous  role  and  is 
most  visible  as  the  slope  of  the  asymptote;  analogous  asymptotic  behavior 
is  insured  by  the  requirement  p > 1.  The  parameter  p,  sensitive  primarily 
it  is  felt  to  deformation,  we  view  as  a shape  factor  controlling  how 
sharply  the  function  rises  towards  the  asymptote**.  The  injection  of 
the  variable  p in  place  of  2 is  inescapably  empirical  but  procures  an 
especially  natural  and  appropriately  conservative  generalization.  Both 
a and  p are  clearly  devoid  of  physical  dimension.  It  is  possible, 
though  not  at  issue  here,  that  for  a large  class  of  situations  the  three 
parameters  can  be  sufficiently  well  prescribed  in  terms  of  relevant 
physical  dimensions  so  as  to  determine  an  explicit  functional  model  of 
predictive  scope  (one  would  for  instance  selectively  borrow  from  rigid 
penetrator  theory  for  a and  vff  formulation)  and  this  theme  may  be  pursued 
later. 


It  is  rather  our  present  purpose  to  advocate  that  form  (2)  serve  as 
a standard  framework  within  which  to  cast  experimental  results;  we 
further  propose  that  explicit  relationships  be  obtained  from  given  sets 
of  (vs>  vr)  data  by  the  systematic  determination  of  triples (a,  v*,  p)  of 
parameters  from  the  regression  procedure  which  is  to  follow  (or,  more 
literally,  from  the  associated  computer  program). 

The  parameter  , marking  the  boundary  of  perforation,  is  often  of 
special  import.  Standardized  adoption  of  the  outlined  scheme  would  in- 
clude acceptance  of  a so  generated  as  being  definitively  the  limit 
velocity  (for  the  given  projectile- target  situation)  implicit  in  the  given 
data  set. 

If 

1.  There  is  resemblance  in  form  to  what  have  been  called  "paver  means," 

p ^ r>  1/p 

expressions  of  the  form  Mp  = (- — ^— ^— ) . ^1/3  for  eszaqple  is  the 

Lorentz  mean  in  the  theory  of  equation  of  state  for  gases. 

It  is  both  interesting  and  irrelevant  to  recall  Fermat’s  notorious 
" last  theorem"  (l?th  Century),  asserting  that  for  a - 1,  vg  > vc  > 0 
and  p > 2,  there  are  no  integral  values  for  the  symbols  vs~,  v ^ and  p 
which  satisfy  (2). 

*AWhile  it  is  of  course  true  that  avs  - vr  0 as  vs  -*•  » for  fixed  p > 1 
( i.e there  is  an  asymptote  of  slope  a emanating  from  the  origin),  it 
is  also  true  that  av  - v •*  0 as  p -*■  * for  fixed  v > v.. 

S X*  5 *. 


13 


4.  Comments,  Perspectives,  and  Implications 


(i)  Initial  expectation  was  that  a non-deforming  projectile  repre- 
sented an  ideal,  that  (mindful  of  traditional  theory)  such  perfection 
could  be  best  accommodated  for  p ■ 2,  and  that  aberrations  should  if 
anything  require  higher  values  of  p (roughly,  that  increase  in  deforma- 
tion should  correspond  to  increased  rate  of  ascent  of  the(vs>  curve 
that  the  "ideal”  should  correspond  to  a limiting  value  of  the  parameter 
p) . Two  seemed  a likely  minimum.  Whatever  the  merits  of  the  reasoning, 
the  expectation  has,  in  our  experience,  been  borne  out  for  situations 
involving  long  rod  penetrators,  but  not  in  general;  fragments  provide  a 
frequent  exception. 

(ii)  As  there  is  occasional  concern  in  the  literature  and  elsewhere 
with  change  or  loss  of  projectile  velocity  effected  by  penetration,  it 
might  be  of  some  value  to  consider  form  (2)  from  this  perspective. 
Letting  A = vg  - vr  we  rewrite  (2)  as: 


v , 0 < v < v„ 

A » 'p  p i/P  • (3) 

K - a <v  - v>  * vs > v£ 

It  is  clear  that  the  corresponding  curve  is  asymptotic  to  the  line 
from  the  origin  of  slope  1 - a;  A is  thus  obliged  to  ultimately,  with 
increasing  striking  velocity,  either  decrease  towards  zero  or  increase 
without  bound  depending  on  whether  a is  one  or  less. 

For  vs  > v and  a < 1 it  is  routinely  verified  that  A achieves  a 
minimum  value  (and  the  corresponding  point  on  the  (v  , v ) curve  has 

p^T  ~1/,p 

slope  1)  at  the  point  v£  = v^  (1  - av  ) . Experimental  affirmation 

that  velocity  loss  can,  as  implied,  decrease  to  a relative  minimum  and 
then  increase  is  noted  e.g.,  by  Goldsmith  and  Finnegan  [1973]. 

(iii)  Another  perspective  can  be  attained  by  considering  transfor- 
mations which  map  the  region  under  the  (vg,  vy)  curve  into  finite  area. 

We  speculate  that  such  device  may  be  of  use  in  providing  additional 
measure  of  parameter  effect.  As  an  example  we  offer  one  such  transfor- 
mation which  maps  functions  of  type  (2)  onto  the  unit  interval  and  yields 
area  as  a function  of  p while  normalizing  the  effect  of  a and  v^. 

v_  v„ 

Let  w = — 1 — and  z = — • For  v > v„  we  rewrite  (2)  as: 
av  v si  v ' 

s s 


w 


(4) 
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The  associated  curve  is  concave,  symmetric  about  the  line  w ■ z and 
contained  by  the  triangle  with  vertices  (0,  1),  (1,  0)  and  (1,  1);  for 
P ■ 2 the  curve  is  a quadrant  of  the  unit  circle.  Area  under  the  curve. 


/‘ 


a 


is  an  increasing  function  of  p with  limiting  values 


of  1/2  and  1 for  1 < p < <*>. 


An  analogous  measure,  and  one  more  computationally  immediate,  is 
provided  by  the  Ly  norm.  Letting  g (z)  - w we  observe  that  g is  an 
element  of  the  (classical  Banach)  space  L (0,  1)*.  The  appropriate 
norm  is  given  by  p 


i 1/p  . 1/p 

*g'  = C/1  lg'P)  = C/  Cl  - zp)  dz) 

P 0 0 


(p+p 


1/p 


Igtp  is  similarly  an  increasing  function  of  p,  and  1/2  < igi^  < 1 for 

1 < p < 00. 


(iv)  This  section  is  concluded  by  noting  an  appealing  property  of 
form  (1)  which  is  not  preserved  in  the  generalization  to  form  (2). 

The  composition  of  two  functions  of  form  (1)  is  itself  of  this  form** 
(composition  is  a binary  operation  on  the  class  of  functions  of  this 
form) , If  (1)  is  regarded  as  the  operable  medium  and  a situation  involves 
a target  comprised  of  several  parallel  plates,  then  (1)  can  be  used 
consistently  for  both  the  whole  target  situation  and  that  of  each 
"subtarget."  In  particular  the  specifying  parameters  for  the  target  can 
be  calculated  from  those  for  the  constituent  plates. 


Unhappily,  such  is  not  the  case  within  the  context  of  form  (2) 
unless  p is  set  at  some  fixed  value.  To  be  quite  specific:  suppose  we 

are  dealing  with  a target  ensemble  of  two  plates  perforated  by  a 


i.e.,  g is  Lebesgue  measurable  on  (0,  1)  and  J*lg  I p It  is  of 

0 

course  equally  true  that  g e (0,  1)  for  all  1 < q < °°  but  in  the 
a..  ! of  p the  observation  is  especially  natural. 

Basically,  because  the  class  of  appropriately  oriented  hyperbolas  is 

V2  2 / 2 2 

x - v and  f „ (;t)  = b Vx  - w 

A z ' 

J~2  2 2 

then  f^  ab  T x - (v  + : the  composition  f^  f£  is  a hyper- 

cl 

bola  of  the  same  type  as  f^  and  f~.  (To  go  further , the  class  of  suck 
hyperbolas  is  an  abelian  group  unaer  the  composition  operation). 


projectile  having  initial,  between-plate  and  final  velocities  denoted 
by  x,  y and  z respectively.  Then  we  would  like  to  say  that 

y * a (x^  - yP)*/P  and  z = b (y^  - w1*)*^  for  some  a,  b,  v,  w,  p and  q. 
But  the  implied  overall  striking- residual  connection  is  then 

z = ab  [if  - vP)^P  - (w/a)^]^^  which  does  not  conform  to  (2).  This 
lack  of  closure  under  composition,  while  not  stri  ly  an  impairment  to 
the  present  concern  of  supplying  effective  representations  of  velocity 
relationships,  does  preclude  full  mathematical  consistency  in  the 
consideration  of  multiple  plate  targets. 


III.  THE  PARAMETERS 


1.  A Regression  Procedure  for  Parameter  Estimation 


We  now  suggest  a mechanism  for  the  determination  of  explicit  form 
(2)  representations  from  given  data  sets.  It  is  convenient  for  the 
moment  to  replace  v , Vf  and  v^  in  (2)  with  x,  y and  c respectively. 


Following  is  a least 

y = 


squares  algorithm  for  fitting  the  form 
0 , 0 < x ^ c 

a (*P  - cV/p,  x > c 


with  constraints  0 a 4 1 and  p > 1 to  a data  set  {(Xj,  y^),  (xn,  yn) ) 

of  pairs  of  non-negative  integers,  each  x.  and  at  least  one  y^  being 

positive.  It  is  assumed  that  n >2,  < xn,  and  x^  < •••  4 xn-* 

Parameters  subjected  to  estimation  are  a,  p and  c (or  otherwise  as 
indicated) . 


For  the  mechanics  of  the  algorithm  we  will  also  suppose  that  p ■<  8; 
this  is  a harmless  restriction  as  there  is  but  minimal  variation  in  the 
form  for  p > 8 (a  and  c fixed)**,  and  it  is  necessary  to  specify  some 
upper  bound  for  p.  Further,  reflecting  upon  the  input  data  and  the 
interval  spacing  selected  for  search  routines,  c will  be  required  to  be 
an  integer  and  p will  be  confined  to  the  set  {1.1,  1.2,  •••,  8}.  These 


*Thai  is,  we  need  at  least  two  data  pairs  and  at  least  two  distinct  x 
values:  additionally , the  data  pairs  are  assumed  to  be  ordered  accord- 
ing to  increasing  values  of  x — this  latter  requirement  will  in  practice 
be  obviated  by  the  inclusion  of  a preliminary  sorting  routine  in  the 
computer  program. 

**A  measure  of  this  variation  is  conveniently  afforded  by  the  L norm 
specified  in  II  4 (iii).  We  note  that  Jgjjg  > .985  while  pl^rn^  J^jj  = !• 
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strictures  have  been  found  reasonable  for  the  particular  application  in 
mind  but  are  of  course  easily  alterable. 

Immediate  interest  is  in  what  will  be  called  the  standard  program, 
SP,  in  which  the  complete  parameter  set  (a,  c,  p)  is  estimated.  There 
is  related  concern  with  the  following  special  options  which  it  is  con- 
venient to  simultaneously  establish; 

0^:  p is  prescribed,  only  a and  c are  determined, 

0,:  a is  prescribed,  only  c and  p are  determined, 

0^:  a and  p are  prescribed,  only  c is  determined. 

In  the  usual  spirit  of  the  least  squares  approach,  we  seek  in  each 
case  appropriately  constrained  parameters  which  approximately  minimize 
the  function 

2 „ „ i/p  2 

* c«  P)  =xf<c  yi  *x2>c  [a  CxiP  ' C?)  ‘ yi]  * 


Root  mean  square  error  (or  standard  error) 


is  expressed  by  4>(a, c,p) . 


With  no  possible  loss  in  precision  (as  measured  by  the  value  of  40 , 
we  require  that  0 < c < and,  for  SP  and  0^,  that  a ~ 1 if  c = xn- 

Oj  and  03 


Let  g(c)  = 2 


x.* 


y. (x.P- 

i i 


cp) 


1/P 


and  h(c)  =2  (x/ 

x.>c 


cp) 


2/P 


c=0, I , . . . ,x 


Then  g > 0 and  h > 0 for  0<c  <x  ; g(x)  =h(x)  = 0;  and 
° n n'  n 


41  (a,c,p)  = 2 y.  - a (2g  - ah). 
i=l  1 


In  the  case  of  0, , if  0 < c < x , ip  is 
1 n r 


a quadratic  in  a having  an 


absolute  minimum  at  a - g/h  > 0 
is  thus  minimized  (exactly)  for 


\ oa 

0 < a < 1 by  taJ 


For  fixed 
min  (g/h. 


c,  ip 


1). 
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4 


I 


Accordingly, 

is 


the  optimal  value  for  a corresponding  to  a given  value  c 

(g(c)/h(c),  g(c)  < h(c) 
a (c)  * j 

g(c)>h(c). 


For  0^  and  0^  respectively,  define: 


(i)  *(c) 


a(c)  [2g(c)  - a(c)h(c)] 


g2  (c)/h(c) , g(c)  < h(c) 


2g(c)  - h(c) , g(c)  > h(c) 


and 

(ii)  <Hc)  =*  a [2g(c)  - ah(c)]. 

n 2 

In  both  cases  the  expression  2 y.  - 4> (c)  needs  to  be  minimized; 

i*l  1 

this  is  accomplished  by  finding  a value  c which  (approximately)  maximizes 
the  function  $(c). 


We  first  generate  an  initial  estimate  cq. 

Let  a = min  {i  = 1,  2,  n : y^  > 0},  xq  * yQ  = 0,  and 


0 = max  {i  = 0,  1,  n;  y^  « 0). 

cq  then  derives  from  a linear  least  squares  fit  of  = ap(xp  - cp)  on 

{(x.p,  y.P)  : i > a}  with  constraints  c > 0,  and  c > x_  if  a > 0.* 
l i o oB 

The  explicit  procedure  for  generating  cQ  is  as  follows: 

Let  r = n - a + 1, 


n 

2 

i=a 


n 

2 

i=o 


x. 

i 


n 

2 

i=a 


n D 

2 (x.  y )p, 

i=a 


In  the  language  of  our  physical  model.  We  are  constraining  the  initial 
limit  velocity  estimate  to  be  at  least  as  large  as  the  largest  striking 
velocity  for  which  there  is  no  perforation  unless  perforation  does  occur 
at  some  smaller  striking  velocity. 

Note  that  the  linear  regression  does  not  use  data  points  to  the  left  of 
the  first  point  for  which  there  is  perforation. 


{ 

\ 


i 
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u * 


and,  to  the  nearest  integer 


0,  rb^  - b2 


lfb4  ' b2  b3 

2 • rbi  ^ b2 


Ibl  - b2 


u > 0 and  b2u  > b_. 


0,  otherwise  . 


Then  let  c 

o 


t,  o f d 

max  (t,  xg),  a > 0 


Normally  c should  be  a very  good  estimate  for  the  parameter  c and, 
as  a simplified  method  suitable  for  desk  calculators,  it  might  be 
acceptable  to  stop  at  this  point;  it  should  then  however  be  realized  that 
c has  not  been  derived  from  a direct  non-linear  least  squares  fit  of  the 
form  to  the  data  but  results  rather  from  a linear  fit  to  a transformation 
of  a subset  of  the  data. 

The  last  step  is  to  search  for  a better  estimate. 

Routine  1:  Starting  from  c = c0,  successively  evaluate  <p(c),  letting 

c vary  by  steps  of  one  unit  in  the  direction  of  increasing  <p  (if  there 
should  be  a local  minimum  at  cQ,  proceed  to  cG+l);  let  c be  the  first 
value  for  c such  that  $(c)  is  at  least  as  large  as  both  Kc-1)  and  ij>(c+l), 
or  such  that  c = 0.  Formally:  let  c = cQ  + jQ,  where 

jo  = 0 if  <Kcq)  > $(co+l)  and  either  cq  = 0 or  4>(co)  > <J>(co-l), 

jc  = min  {j  > 1 : $(c.o+j+l)  < <t»(co+j)}  if  <Kcq)  < 4>(co+l), 

jQ  = -min  (j  > 1 : 4> Ccq' j -1)  < <Kco-j)  or  j = cq}  if 

♦ CV«  * *(c0)  < ♦(c0-l). 

While  there  is  no  assurance  that  c provides  the  largest  possible 
value  for  d> Cc) , it  does  better  than  the  initial  estimate,  cQ,  which  is 
itself  expected  to  be  good  and  further  searching  is  considered 
unwarranted . 

For  the  parameter  estimates  are  a = a (c)  and  c = c. 

For  the  estimated  parameter  is  c = c. 
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SP  and  0, 


We  rely  upon  much  of  the  preceding  technique  while  allowing  p to 
vary. 

^ f 

For  a fixed  value  of  p,  let  g(c),  h(c)  and  a( c)  be  as  before  and  let 
$(c)  be  defined  as  in  (i)  and  in  (ii)  for  the  cases  SP  and  Qj  respec- 
tively. 

For  p * 2,  determine  cc  as  described  previously;  then  use  Routine  1 
to  find  the  corresponding  c,  which  we  now  denote  as  c(2). 

For  p = 2.1,  2. respectively,  let  c(p)  be  obtained  from  its 
predecessor  by  using  Routine  1 with  cQ  replaced  by  c(p  - .1).  Similarly, 
for  p « 1.9,  1.8,  •••1.1,  let  c(p)  be  successively  generated  by  using 
Routine  1 with  cQ  replaced  by  c(p  .1). 

Let  $(p)  * $(c(p)),  P - 1.1  1-2,  ”-,8. 

Approximate  maximisation  of  $(p)  is  desired  and  we  follow  a search 
routine,  starting  from  p = 2.*  Explicitly:  let  p = 2 ♦ k/10  where 
k * min  {j  = 0,1, •••,60  : $(2+j’10)  > $(2.1+j/10)  or  j » 60}  if 
$(2)  < $(2.1); 

k = - min  (j  = 0, 1.  • • • ,9:$(2-.i/10)  > $(1.9-j/10)  or  j = 9} 
if  $(2)  > $(2.7). 

For  SP  the  parameter  estimates  are  a = <x(c(p)),  c * c(p)  and  p = p. 

For  0^  the  estimates  are  c = c(p)  and  p = p. 

2 . A Remark 


Reliance  upon  a direct  non-lines.r  regression  technique  necessitates 
high-speed  computing  ability.  The  fora  does  not  seem  to  lend  itself  to 
linearization  in  any  acceptable  manner  but  we  do  at  least  expect  to 


Starting  from  p = 2 is  primarily  a matter  of  computational  expedience; 
we  expect  that  in  many  instances  the  optimal  value  of  p will  be  near  two. 

The  search  works  away  from  two  in  the  direction  of  increasing  $ and  stops 
when  a local  maximum  is  detected.  There  is  no  general  guarantee  of  hav- 
ing hit  upon  the  absolute  maximum;  however  we  have  tried  more  extensive 
routines  on  a large  and  diverse  class  of  data  sets  and  have  by  comparison 
for'nd  the  simple  (and  efficient)  routine  given  here  to  be  consistently 
satisfactory  (typically  the  function  $,  for  the  given  data,  had  but  one 
maximum. ) We  are  confident  i.iat  this  routine  does,  for  the  intended 
application , provide  a good  estimate. 
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benefit  in  terms  of  accuracy  from  a direct  approach.  Prevailing  against 
linearization  are  the  nonlinearity  of  the  form  in  the  parameter  p and 
the  special  role  of  the  parameter  (or  c) . 

If  one  were  content  to  take  form  (1)  (or  (2)  with  p prescribed), 
restrict  this  form  to  its  support,  delete  certain  data  points  and  trans- 
form the  remainder,  then  linearization  is  feasible.  A problem  here, 
besides  loss  of  generality  and  incomplete  utilization  of  information, 
is  that  a sound  decision  on  point  deletion  can  be  critical  and  is  not 
always  apparent.  An  obvious  but  not  consistently  judicious  choice  is 
to  delete  all  points  for  which  there  is  no  perforation;  i.e.,  the  linear 

2 2 2 2 2 2 
formvr  = a (v3  - vf  ) is  adapted  to  those  points  (vs  , vr  ) for  which 

vr  is  positive.  Such  procedure  is  recurrent  in  the  literature  and  goes 

back  at  least  to  Robertson  [1943].  A rather  more  sophisticated  approach 

along  these  lines  was  used  for  generating  initial  estimates  in  the 

preceding  algorithm. 

For  the  price  of  nonlinearity  and  a higher  level  of  machine  depen- 
dence we  find  it  far  preferable  to  systematically  fit  a complete  and 
more  general  form  to  a complete  data  set. 

3.  Some  Statistical  Contingencies 

The  effort  here  is  to  suggest  some  plausible  areas  for  statistical 
exploration.  The  questions  raised  will  likely  require  better  formula- 
tion from  keener  perspectives  at  another  time  and  we  will  not  at  this 
point  be  particularly  inhibited  by  lack  of  precision. 

Till  now  we  have  adopted  a deterministic  view  of  velocity  dependence. 
We  have  assumed  the  pretense  of  being  able  to  deal  with  a specific 
residual  velocity  consequent  to  a given  striking  velocity;  more  real- 
istically we  should  perhaps  enquire  of  the  probability  that  for  a given 
striking  velocity  the  residual  velocity  (regarded  as  a random  variable) 
will  belong  to  a specified  interval.  It  is  transparent  that  for  each 
vs,  Prob  [0  < vr  < vs]  = 1;*  significant  statements  of  this  nature  would 
however  be  welcome. 

To  what  extent  can  meaningful  assessments  of  confidence  be  attributed 
to  estimated  parameters?  Size  of  the  data  set  is  clearly  important. 

Also  significant  surely,  and  in  different  ways  for  the  different  parame- 
ters, is  the  distribution  of  striking  velocities;  a concentration  of 
p unts  about  the  estimate  v^  is  bound  to  enhance  confidence  in  that 
estimate  but  reflects  little  on  the  parameter  a which  is  mostly  influ- 
enced by  points  with  high  striking  velocity.  Estimation  and  exploitation 
of  relevant  variance  measures  is  probably  essential. 


Unless  fortuity  in  measurement  is  considered,  to  admit  of  the  possibility 

that  v > v . 
r s 
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The  parameter  v£  is  often  of  primary  interest.  One  would  ideally 
hope  that  some  reasonable  assumption  about  probability  distribution, 
interacting  with  statistical  invention,  could  be  brought  to  assert  that 
for  a reasonable  data  see  the  quantities  Prob  [vy  = 0 j vs  < - e]  and 

Prob  [Vj.  >0  | vs  > + e]*  are  reasonably  large  (proximate  to  one)  for 

reasonably  small  positive  values  of  c.  It  would  be  additionally  provi- 
dent were  the  above  probabilities  to  equal  one-half  for  e «*  0,  in  which 
case  vK  would  be  equivalent  to  what  has  elsewhere  been  termed  "V50",  the 
concept  of  which,  along  with  other  piobabilistic  phenomena,  has  provoked 
the  perpetration  of  much  nonsense  in  ballistic  literature. 

Of  allied  interest  would  be  a measure  of  relative  confidence  in 
different  v^'s  extracted  from  different  data  sets.  Consider  the  follow- 
ing situation.  A collection  of  v^'s  is  at  hand  and  one  is  equipped  with 
a general  form  purportedly  able,  over  some  broad  range  of  situations 
including  those  which  yielded  the  v^'s,  to  represent  limit  velocity 
in  terms  of  the  physical  set-up  (materials,  geometry,  etc.);  e.g., 
variants  of  the  de  Marre  form*  have  been  fashionable  in  this  regard  for 
nearly  a century.  One  then  regresses  to  the  collected  v^'s  to  evaluate 
parameters  - but  ideally  the  various  v£'s  should  be  variously  weighted 
in  the  relevant  regression  function  so  as  to  reflect  their  various 
relative  degrees  of  reliability. 

Another  problem  involves  the  design  of  methodology  for  data  acquisi- 
tion relative  to  a given  projectile-target  situation.  With  reference 
to  a priori  value  judgments  about  information  desired,  and  constrained 
by  economic  and  physical  limits,  a suitable  and  efficient  experimental 
strategy  needs  to  be  ordained.  A new  dimension  of  complexity  is  injected 
at  this  stage  by  the  introduction  of  another  random  variable;  striking 
velocities  can,  we  suppose,  be  regarded  as  deterministically  Treasured, 
but  they  certainly  cannot  be  so  controlled. 


* . • * 

Fy'ob  [A|B]  is  the  (conditional) probability  measure  of  A gzven  B. 

The  most  prevalent  de  Marre- type  form  and  one  in  current  local  usage , 
a dimensionally  purified  revision  of  the  original , can  be  written  as 

v^  = (t/d)Y  yj kd’Vm  where  d and  m are  projectile  diameter  and  mass , and 
t is  target  thickness  (or  thickness  times  a function , such  as  secant , 
of  incidence  angle),  k,  with  units  of  force  per  area  (e.g.  a multiple 
of  yield  stress) , and  y are  the  parameters  to  be  estimated.  A 
(linearized)  regression  procedure  in  this  case  is  apparent , The  form 
is  (perhaps  obviously)  empirical  but  not  without  an  element  of  physical 
appeal , especially  in  the  rigid-projectile  context  for  which  it  was 
originally  intended ; indeed  some  quite  serious  theoretical  analyses, 
invoking  principles  of  elasticity  and  hydrodynamics,  have  produced 
models  of  this  type  with  values  for  y of  1/2  and  1 (our  data-dependent 
least  squares  determinations  for  y have  typically  been  near  .8). 
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If  liait  velocity  is  the  sole  objective,  the  Attempt  should  be 
towards  obtaining  perforations  and  non-perf orations  within  a snail  range 
of  striking  velocities.*  Should  the  intent  be  for  uniforaity  of  confi- 
dence in  the  full  relationship,  or  if  there  is  a weighted  concern  with 
perhaps  disproportionate  interest  in  liait  velocity,  then  an  iaposing 
variety  of  concepts  needs  to  be  formalized  and  an  appropriate  experiment 
designed.  At  issue,  briefly,  is  the  prescription  of  standard  procedures 
for  generating  sets  of  (vs,  vp  data  from  which  liait  velocities  and 

overall  velocity  correlation  can,  with  measurable  reliability,  be 
systematically  and  definitively  determined. 


IV.  THE  COMPUTER  PROGRAM 

1.  Introduction  and  Examples  of  Output 

The  algorithm  of  Section  III.  1.  is  encoded  in  a Fortran  Program 
to  generate  parameter  estimates  from  a given  data  set.  Plot  routines 
provide  separate  graphic  displays  of  each  of  the  consequent  form  (2) , 
(3)  and  (4)  representations;  the  form  (2)  representation  (the  (vs>  v ) 

curve)  being  no  doubt  of  predominate  interest.  In  each  case  the  given 

data  set  or  the  appropriate  transformation  thereof  is  also  graphed;  for 

(2)  and  (3)  relevant  asymptotes  are  plotted  as  dotted  lines.  In 

deference  to  machine  notational  limitations,  the  symbols  v . v , v,,  a 

s r <c 

and  p are  replaced  by  X,  Y,  C,  A and  P respectively  in  plots  and  print- 
out. S will  designate  the  root  mean  square  error  (c.f.  page  17) 
associated  with  the  fit  of  the  determined  (v  , v ) curve  to  the  given 

(v.,  v ) data  set. 
s r 

Although  our  special  concern  in  this  report  is  with  the  standard 
program  SP  (in  which  all  three  parameters  are  determined),  we  provide 
as  well  for  the  options  0^,  02,  and  0^  described  earlier  and  addi- 
tionally for  the  possibility  of  plotting  graphs  and  relevant  points 
when  all  three  parameters  are  initially  prescribed. 

One  page  of  computer  printout  includes  a tabular  listing  of  the 
data  set  (observed  striking  and  residual  velocities)  along  with 
some  derived  quantities  which  may  be  of  peripheral  interest  (e.g., 

*One  eonjeolwea  that  a so  generated  is  inherently  "better"  than  the 

analogous  extraction  from  Sensitivity  Analysis,  in  which  there  is 
sensitivity  only  to  quantal  response  (yes  or  no  as  regards  perforation ) 
and  which,  provided  residual  velocities  are  watched,  seems  far  less 
exploitative  of  information. 
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the  "error"  at  each  observed  striking  velocity  between  the  observed 
residual  velocity  and  the  residual  velocity  associated  with  the  determined 
fit) . 

For  SF  and  0^,  a second  page  of  printout  tabulates  "optimal" 

values  for  A and  C corresponding  to  each  value  for  P considered  by 
the  search  routine  described  on  page  20;  the  root  mean  square  error  S 
associated  with  each  such  P is  also  calculated  and  printed.  In 
particular,  the  "optimal"  form  (i)  fit  is  determined  (recall  that 
the  search  routine  in  question  starts  from,  and  hence  always  examines, 
the  case  P *2). 

In  the  following  pages  we  offer  examples  of  machine  output  as  generated 
by  the  standard  program.  These  are  meant  to  be  more  illustrative  than 
substantive;  actual  projectile-target  situations  involved,  though 
partially  and  cryptically  specified  in  titles,  are  of  little  concern  and 
we  do  not  elaborate. 

Figure  1 provides  displays  of  striking  versus  residual  velocity 
(form  (2)  representations)  generated  by  two  different  data  sets.  Figure 
2 depicts  striking  velocity  versus  residual  velocity  and  striking  velocity 
versus  velocity  loss  (form  (2)  and  (3)  representations  respectively) 
associated  with  the  same  data  set.  Next  is  an  example  of  the  complete 
output  (in  actual  output  size)  from  a single  set  of  input  data:  tables 

la  and  lb  comprise  the  computer  printout;  figures  3a,  3b  and  3c  are 
respectively  the  forms  (2),  (3)  and  (4)  plotter  displays. 
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Figure  1.  Two  Examples  of  Generated  vg  versus  Curves  (Form  (2) 
Representations) 


Table  la.  Sample  of  Computer  Printout:  First  Page 
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AND  Y*  ARE  IN  METERS/SEC 


Table  lb.  Sample  of  Computer  Printout:  Secord  Fage 
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Sample  of  Plotter  Output:  v_  versus 
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Figure  3b.  Sample  of  Plotter  Output: 


2.  Input  and  Program 

In  Appendix  A is  found  a guide  to  the  input  for  running  the  regres- 
sion procedure  program  listed  in  Appendix  B.  The  program  is  written  in 
Fortran  language  for  BRLESC  but  could  be  adapted  to  other  computers , 
with  a few  minor  changes,  i.e.,  replacements  for  SORTXY,  SCOOP,  length  of 
Hollerith  statements  and  USE  (MAIN  C(**40N).  SORTXY  (X,  Y,  N)  sorts  the 
elements  of  vector  X into  non- decreasing  order  and  moves  the  elements  of 
Y so  they  correspond  to  the  original  values  of  X.  SCOOP  is  the  BRLESC 
implementation  of  the  basic  Fortran  software  package  described  in  the 
publication  "Programming  CALCOMP  Pen  Plotters"  [publication  No.  1006A, 
California  Computer  Products,  Anaheim,  California,  1969].  A word  of 
caution:  XL,  the  third  argument  in  PLOTS,  is  the  length  of  paper  required 
instead  of  the  logical  output  device  number  as  commonly  used.  This 
program  is  set  up  for  30  inch  width  paper. 

Any  or  all  the  options  and  alternatives  can  be  employed  by  setting 
the  appropriate  flag  for  each  data  set.  The  plots  and  tabulation  will 
indicate  which  parameters,  if  any,  are  specified.  In  addition  the  first 
and  third  plots  can  have  different  scales.  Hie  usual  procedure  is  for 
SCOOP  to  determine  the  initial  X and  Y values  and  data  units/inch  but 
this  can  be  bypassed  with  the  flag  FIXAX  and  the  alternative  values 
stipulated. 

Hie  program  employs  13K  memory  and  a collection  of  sample  problems 
averaging  fifteen  pairs  of  data  points  required  an  average  of  19  seconds 
of  computer  time  per  problem. 
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APPENDIX  A 


INPUT  DATA 


AFRWMS  A 


INPUT  DATA 


CARD 

TOftMAT 

CO  LOWS 

PROGRAM  NAME 

UNITS 

REMARKS 

i 

FIO.0 

I- 10 

XL 

inchet 

Maximum  amount  of  plotting  pepor  for 
this  computer  run.  (Approximately  IS 

in  chat  per  set  of  3 plots) 

F10.0 

11-30. 

PINCH 

Value  used  for  incrementing  parameter  F, 
N.l  is  used  in  this  report 

2 

7AIO 

1-70 

TITLE  (I) 

Title  to  head  tables  and  graphs  Cup  to 
70  extractors) 

3 

IS 

1-S 

N 

Number  of  data  points.  Dimension  state- 
ments are  set  up  for  97  hut  can  be 
altered 

IS 

6-10 

IFF 

Dimension  flag  for  (X,  Y)  data  (m/sec*0, 
ft/*ec«l) 

IS 

1 1- IS 

FIXAX 

Axis  flag  for  Plot  1 and  Plot  3 (not 
specified«0.  specif ied“ 1) 

IS 

16-20 

PLT 

Plots  only,  where  P,  A and  C are  given 

IS 

21-2S 

OPT(l) 

Flag  for  parameter  P given  (ycs-1, 
no-0) 

IS 

26-30 

0PT(2) 

Flag  for  parameter  A given  (yes-1, 
no-0) 

IS 

3I-JS 

OPTO) 

Flag  for  parameters  P and  A given  (yos«l 
no-0) 

IS 

36-40 

0PT(4) 

Flag  for  eliminating  standard  program 
(yes-1,  no-O) 

4 ate. 

•F10.0 

1-M 

X(I).  Y(l) 

ft/tcc  or 

(vs,  vr)  data,  four  pairs  to  a card,  in 

■/sec 

units  of  m/sec  or  ft/sec  as  indicated 
by  IFF  flog 

HU;  FOLLOWING  CAROS  ARC 

DLPENOENT  ON  ABOVE  FLAGS: 

IF  NONE, 

IGNORE 

If  FIXAX-1 

F10.0 

1-10 

PLFX1 

■/sec 

Minimum  X value  for  Plot  1 horizontal 

axis 

F10.0 

11-20 

PUlXl 

■/sec 

Number  of  data  units  per  inch  for  Plot  1 
horizontal  axis 

F10.0 

21-30 

PLDY1 

a/sec 

Number  of  data  units  per  inch  for  Plot  1 
vertical  axis  which  starts  at  zero 

F10.0 

31-40 

PLOXJ 

■/see 

Number  of  data  units  per  inch  for  Plot  3 
horizontal  axis  which  starts  at  zero 

F10.0 

41-50 

PL0Y3 

■/sec 

Number  of  data  units  per  inch  for  Plot  3 
vertical  axis  which  starts  at  zero 

If  PLT-1 

F10.0 

1-10 

C 

■/sec 

Parameter  of  equation 

F10.0 

11-20 

A 

Parameter  of  equation 

FJO.O 

21-30 

P 

Parameter  of  equation 

If  OPT(I)«i 

F10.0 

1-10 

P 

Parameter  P given  J One  card  per  option 

I-l.  3 

| requested.  For  on 

FlOaO 

11-20 

A 

Parameter  A given  j (1)  P only,  OPT  (2) 

l A only,  and  OPT  (3) 
1 both  P «nd  A 
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oftoflo  oononnot'ioono  o norm 


dimension  osaoo),  ovuooj,  eyuoo),  optia),  yciiooi 

COMMON  k,  AX ( 100  I t AY ( 100 ) f BDXI5S,  BOY ( 5)  , 8DZ(5),  BX<7) , (XY(7t, 

1 CX  ( 103 ) . CY(I3C),  FIXAX,  FX(IOO),  L»  N,  P,  PIHCR,  PtNV,  S,  0 VStJ  • 

2 TA(  3)  • TB  ( 3 ! * TC(2),  TD(2),  TFC2),  T0(3),  TITLEITJ,  XIIOO). 

3 Y( 107 ) , PLFX1,  PL0X1,  PLFYl,  PLDYl,  PLFX3,  PLDX3,  PLFY3.  PL0Y3, 

4 PLT , IUUF(  10001,  XINCj  YINC,  TE,  C,  TGI2I,  TH(2),  NCASE 

INTEGER  ALPHA, BETA, F I X AX , OP T * PL T 

DATA  (BDXtn.BDYUJ.BOZm,  1-1,51  /-l. 2,  -1.5,  -1.8,  9.3,  -1,5,  -i 
1*8 » 9*3,  6.5 , 6.2,  — 1.2,  6.5,  6,2,  —1.2,  —1,5,  —1.8/ 

DATA  (Bxm.bYU), 1 = 1,7)  / 5. 87 , 2.46,  5.76,  2.46,  5.76,  2.25,  5.69 
1,  2.215,  5.76,  2.18,  5.76,  1.97,  5.87,  1.97/ 

OATA  (TA( II,  1 = 1,31  / 10HX  STRIKING,  1 OH  VELOCITY  , lOHCM/SECI  / 

OATA  (Tam, 1-1,31  / 1CHY/AX  = U-(C,  ICH/X)  5 , X>,  IGHC  / 


OATA  ITCIII,M.2I  /19HA(X  -C  I , 19H,X>C  / 

OAT  A (Tom,  1*1,21  / 10HP  P , 1CHI/P  / 

OATA  TE  /7H1GIVENI/ 

OATA  (TF( I), 1*1,2)  /10HX-AIX  -C  ),  10H  ,X>C  / 

OATA  |TC(1). 1*1,2)  / 10HX-Y  (M/SEC,  10H)  / 

OATA  ITOm,  1*1,3)  / 10HY  RESIOUAL,  10M  VELOCITY  , 1CHCM/SEC)  / 

OATA  (THI I), 1*1,2)  /1QHPARAMETEKS,  10M  GIVEN  / 

READ  (5,37)  XL,PINCR 


XL  - TOTAL  LENGTH  OF  PAFER  NEEDED  FOR  GRAPHS 

PINCR  - INCREMENT  FOR  P FOR  OPTIONS  2 AND  4 

CALL  PLOTS  ( IBUF, 1000, XL) 

CALL  PLOT  10,0,-36.0,-31 
CALL  PLOT  (5. 0,2. 0,-3) 

X INC* 15. 

VI NC-9. 

FKIll'O. 

OO  1 1*2,50 

1 FXm*FLOAT{  I-U+.02 
FX ( 5 1 ) *1,0 

FX ( 52 ) =0.0 
FX(5 j)=»25 

*♦**  TITLE  WILL  HEAD  ALL  GRAPHS  AMD  TABULATIONS 

2 READ  (5,35)  ( T I TL t ( I ) , I* 1 , 71 
WRITE  16,391  (TITLE! I) ,1*1, T) 

IHEAOG-Q 

READ  (5,36)  N, IFF, FIXAX, PLT, (CPT(I), 1 = 1,4), (XU), Y(J), J-l.NI 

N - NUMBER  OF  DATA  PAIRS  ( D l HEMS  I ONI NG  ALLOWS  MAXIMUM  OF 

IFF  - DIMENSION  FLAG  FTR  OATA  M/SEC=0  FT/SEC-l 

F IX  AX  - FIXED  SCALE  FOR  AXIS  *<0*0  YES*1 

PLT  - PLOT  WITH  GIVEN  A,P,ANO  C YES-1  N0=0 

OPTU)  - P GIVEN  YES-1  N0=0 

OPT? 2)  - A GIVEN  YES-l  H0=0 

uPTij)  - A AND  P GIVEN  YES*I  N0=0 

OPT ( 4 1 - ELIMINATE  STANDARD  PROGRAM  YES-1  NO-O 

X - STRtKlNG  VELOCITY 

y _ resiouAL  VELOCITY 

IF  ( FIXAX. GT .0)  READ  (5,37)  PLFX1,PLDX1 , PLDYl, PLDX3,PLDY3 

PLFXl  - STARTING  VALUE  ON  HORIZONTAL  AXIS  PLOTl 

PL OX l - DATA  UNITS  PER  INCH  HORIZONTAL  AXIS  PLOTl 

PLDYl  - OATA  UNITS  PER  INCH  VERTICAL  AXIS  PLOTl 

PL0X3  - OATA  UNITS  PER  INCH  HORIZONTAL  AXIS  PL0T3 

Preceding  page  blank 
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OOUOU  O -J  o o 


C PLDY3  - DATA  UNITS  PER  INCH  VERTICAL  AXIS  PLOTS 

IP  UFf.EQ.QI  GO  TO  4 
C **•*  CONVERTS  FEET  TO  HETERS 

00  3 K- 1 ,N 
X(K1«X(K£*,3048 

3 r<KI-Y(Kl*.304a 

C »♦*  ARRANGING  OATA  IN  INCREASING  ORDER  FOR  STRIKING  VELOCITIES 

4 CALL  SORT (X&Y,NI 
SYSQ-O. 

DO  5 I-l.N 

5 SYSQ«SYStt*YI  H*YI  II 
IF  IPLT.GT.06  p.cAD  C S* 37)  C,A,P 

C - LIHIT  VELOCITY 

A - SLOPE 

P - EXPONENT  OF  FIT 

WRITE  (6,331  (Xm,Y<I)»I«l.NI 

6 DO  32  1*1*4 
IF  (PLT.GT.Oi  GO  TO  25 
IF  ( CPT  ( H. EQ.O.  AND. L.LT.4 1 GO  TO  32 
IF  ( L.EQ. 4.  4ND.0PT ( 4 ) *F.Q«  1 ) GO  TO  32 
IF  (L.EQ.4)  GO  TO  8 
JFIN-C 
JOPT-1 

7 IF  (OPT(L).LT.JOPT)  GO  TO  32 
READ  (5,37)  P,A 

P , A - PRESCRIBED  VALUES.  NONE  TO  THREE  DATA  CARDS  DEPEND- 
ING ON  OPTIONS 

JOPT-JOPT+l 
GO  TO  9 

8 JF IN- I 

9 00  10  I-UN 
ALPHA- I 

IF  (Y(I).GT.O.)  GO  TO  II 

10  CONTINUE 

11  UETA-0 
NN1-N-1 
CO  12  I-UNfll 
IF  (Y(I).EQ.O.)  8ETA-I 

12  CONTINUE 

C **♦*  INITIAL  ESTIMATE  OF  C 

81-0. 

B2-0. 

83- 0. 

84- 0. 

IF  (L.EQ.2.0R.L.EQ.4)  P-2. 

DO  13  K>1»N 

IF  tYUl.GT.O.)  GO  TO  14 

13  CONTINUE 

14  03  IS  I -K *N 
XP-XI I»#*P 
YP-Y( I 1**P 
Bl"BL«XP**2 
82-82+XP 
B3-B3*YP 

15  84“B4*XP*YP 


42 


R«.N-K*1 
Rft l»R*Bl 
82SQ«B2*B2 

IF  (RB1.EQ.B2SQ)  CO  TO  16 
U«(R*B4-02*B3»/IR&:-B2SQI 
GO  TO  17 

16  U*0. 

17  IF  (U.LE.O.  i GO  TO  IB 
T».(B2*U-B3)/(R*UI 

IF  (T.LE.O.)  GO  TO  18 
T«T**(  l./P) 

GO  TO  19 

18  T»0. 

19  TI 3 A INTI T ) 

TK-T-TI 

T-TI 

IF  ( TK«GT » « 5 ) T«TI*l. 

IF  ( ALPHA. L c. BET A ) CZ* AM AX 1( T,0. I 

IF  ( ALPHA. GT. BETA)  C Z* AMAX 1 ( T , X < BE TA )> 

C-C  l 

WRITE  (6, AO) 

Il-O 

IF  I IHEACC.eC.l)  GO  TO  34 
IHEADG=l 

20  CALL  PARAH6TERS 
WRITE  <6.91)  P»A,C,S 

IF  (L.EQ.2.0R.L.EQ.4)  GO  TO  26 

21  WRITE  (6,42) 

C 4#**  SUMMARY  TABLE  OF  DATA  FOR  PLOTS 

AVG»0. 

PINV’l./P 
DO  24  I=WN 

IF  { X ( D.LE.C)  GO  TO  23 
YC ( I ) 3 A* ( X I I )**P-C**P)**P5NV 
-2  EY(I)"Y(  I).-YC(  I) 

DS(  I )»X(I  l-YC( I ) 

0Y(1)=X( I)-Y(I) 

AVG3AVC*EY(I) 

GO  TO  24 

23  YC(I>*0. 

GO  TO  22 

24  CONTINUE 

AVG* AVG/ FLO  AT ( N ) 

WRITE  (6,43)  (TITLE( 1 1 , 1 -1, 7) ,P, A,C, AVC 

WRITE  (6,44)  (XI  I ),Y( I ),YC( I ),EY( I ),DY(I) ,DS(I ) #1*1, NJ 

WRITE  (6,45) 

IF  (L.EQ.ll  WRITE  (6,46) 

IF  (L.EQ.2)  WRITE  (6,471 
IF  (L.EQ.3)  WRITE  (6,48) 

WRITE  (6,42) 

25  CONTINUE 

IF  (PLT.GT.O)  L»0 

C 4*4*  PLOT  * VS  Y (STRIKING  VELOCITY  VS  RESIDUAL  VELGCITY) 

CALL  PLOT l 

C •***  PLOT  C/X  VERSUS  Y/AX 

CALL  PLOT 2 

C 444*  PLOT  X VS  X-Y 

CALL  PL0T3 

IF  (PLT.GT.O)  GO  TC  33 
IF  IJFIN.EQ.O)  GO  TO  7 
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o r.  o 


C 


C 


CO  TO  32 

♦ *«  QETERN INAT  ION  OF  P 

26  IF  (II. or. 01  CO  TO  28 
1 1 -1 

0IR«PINCR 

KN*0 

27  PREVA-A 
PR£VC*C 
PREVS-S 
PREVp.p 
P*P*0!R 
CO  TO  20 

28  IF  (S.Ot.PREVSI  GO  TO  30 

IF  4P.CE.3.0.0R.P.LE.1.1)  GO  TO  21 

29  MM-l 

GO  TO  27 

30  IF  (MH.GT.O)  CO  TO  31 
OlR*-PlNCR 

GO  TO  29 

31  C-PREVC 
A*PREVA 
P-PREVP 
S-PREVS 
CO  TO  21 

32  CONTINUE 
GO  TO  2 

33  Pt  T *0 
IH£ ADC* 1 
CO  TO  6 

34  WRITE  (6,39)  (TITLE!  I), 1*1,7) 

WRITE  (6,38)  <X(I),Y(l),I*l,N) 

GO  TO  20 


35  FORMAT  (7A10) 

36  FORMAT  4 d 1 5/ ( 8F 10.0) ) 

37  FORMAT  (BF10.0) 

38  FORMAT  (//6X, 1HX.9X, 1HY.9X, IHX , 9X, IHY, 9? - IHX, 9X , 1HY.9X, 1HX „9X» 1HY/ 
l/f3X«F6.Qr4X,F6«0,3X»F6.0,4X,Ffc»C»4X,F6.0,4X,F6»Q»4X,F6»0,4X»F6.0) 


2) 

39  FORMAT 

40  FORMAT 

41  FORMAT 

42  FORMAT 


(1X.7A10) 

(///) 

4 ' P* 1 , F5.2 , 3X , ' A 
« IH1) 


' ,F6.3,3X,*C»*,F6.0,3X,'S-,,F9.3) 


■» A rw«\r  ( inti 

43  FORMAT  4//,lCX,7A10////l5X,2HP*tF4.l,5H»  A*,F5.2,5H,  C*»F 

1 M/SEC,  AVG(Y-Y,1*,F6.0/15X,49HY**A4X**P-C**P>*'M1/P)  IF  X> 
20  OTHERWISE  /// 17X , IHX , 8X , IHY, 6X, 2HY * ,6X .4HY-Y® .5X.3HX-Y, 
3V  ) 


,F6.0,16H 
X>C,  Y*» 
6X, 4HX- 


JY'  ) 

44  FORMAT  4 14X , F6.0 , 3X , F6.0 , 3X , F6.0, 3X, F6, 0, 3X , F6.0 ,3X,F6.G> 

45  FORMAT  4 /// ICX , 28HX • Y AND  Y'  ARF  IN  METERS/SEC,//) 

46  FORMAT  4 10X, •PGIVEN* ) 


47 

48 


FORMAT  4 LOX , • AGIVEN* ) 
FORMAT  4 10X, 'PANDAGIVEN' ) 


SUBROUTINE  PARAMETERS 
COMMON  (USE  MAIN) 

DIMENSION  XP4100),  YP4100) 

M-0 
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onn 


C «***  DETERMINATION  OF  A 

00  1 K»liN 
XP (K) *XI Kl **P 

1 YPtlU«YIKl*«'P 
PINV-l./P 

2 IF  (C.LT.O.)  CO  TO  0 
CP“C**P 

00  3 K»1,N 

IF  (X(KI*G£«C)  GO  TO  A 

3 CONTINUE 
A G*.Q. 

H«.U. 

DO  5 I«K.N 

XPMCP»(XPIII-CP  J *YP INV 
6«G>Y ( 1 1 *XPHCP 

5 H»H+XPHCP**2 
GDH*G/H 

IF  IL.EO.l.OR.L.EO.A)  A-AMINHGDH,  1.  ) 
PMI=A*(2»*G-A*H) 

C **  DETERMINATION  OF  C AND  S 

IF  ( M.GT . 0 ) GO  TO  7 
H*1 

OlR«l* 

MM*0 

6 SAVE  A" A 
SAVEP-PHI 
SAVEC*C 
SAVEG»G 
SAVCH»H 
C*C»OIR 
GO  TO  2 

7 IF  (PHI.GE.SAVEP)  GO  TO  8 
IF  (MM.GT.OI  GO  TO  9 
OtR—  l. 

8 HM-l 

GO  TO  6 

9 A-SAVEA 

C-SAVEC  , _ 

PS' l3  AOS  I SY$0-A*<  2.*SAVEG-A*SAVCH>  ) 
S*SCRT(PSI/FLOAT(N) ) 

RETURN 

C 

END 


SUBROUTINE  PLOT! 

COMMON  (USE  MAIN) 

PINV-l./P 

IF  (FIXAX.GT.O)  GO  TO  16 
IF  (C.LT.XI l ) ) GO  TO  13 

• fdi  cm  C I V.R.O*Nt  1) 

1 ASSURES^ Y-AXIS  STARTS  AT  ZERO 
Yl N*l 1*0. 

N«N+1 

CALL  SCALE  IY,5.0»N,1I 
YI NS*VIN*1) 

Y(N*1I-TIN*2I 

N«N-1 
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IF  IJJ.EQ.l)  GO  TO  14 

2 AA.C-X(N*2) 

IF  IAA.LT.X(N+lt)  GO  TO  15 
C ***  PLOT  TITLES. AXES,  AND  PARAMETERS 

3 CALL  AXIS  (0.0,0.0,TA,-30,8.0,Q.CSX(NUI,X(N*2)) 

CALL  AXIS  <0.0,0.0,TO,30,5.0,90.0,Y(N*l) ,V(N*2)> 

IF  (FlXAX.GT.O)  GO  TO  19 

CAU  LINE  (X.Y.N.l.-ltl) 

* CALL  SYMBOL  (-.5, 5. 5, *14, TITLE, O.O, 701 
CALL  SYMBOL  ( 6.46, 2*32, .07, 5H<  <t 0.0,5) 

CALL  SYMBOL  ( 5.9, 2. 25, . 14, 8H0,  0 X C.0.0,8) 

CALL  SYMBOL  I 6.46, 2. 25, .07, 5H*  >,0*0,5) 

CALL  SYMBOL  ( 5.34, 2. IB, . 14,2hY*,0.0,2) 

CALL  SYMBOL  ( 6. 32, 2. 145, .07, TD,G.O, 20) 

CALL  SYMBOL  I 5, 9 , 2.04, . 14 , TC ,C. 0, 20 ) 

CALL  SYMBOL  ( 5. 34,  1. 59, . 14, 2HP= , 0.  0, 2 ) 

CALL  NUMBER  ( 5. 62, 1 . 59 , . 14, P ,C. 0,1 ) 

IF  (L.EQ.  1.0R.L.CQ.3)  CALL  SYMBOL  (6.32,1. 59,.  14, TE, 0.0, 7) 
CALL  SYMBOL  ( 5. 34, 1. 3C , . 14, 2HA*,0. 0,2) 

CALL  NUMBER  ( 5. 62 , 1. 36 , . 14, A.C.O, 2 ) 

IF  (L«  EQ. 2.0R.L . HQ. 3 I CALL  SYMBOL  (6.32,1.38,. 14, TE, 0.0, 7) 
CALL  SYMBOL  ( 5. 34,  1. 17,  . 14,  2HO,  0.  0, 2) 

CALL  NUMBER  (5.62, 1. 17, . 14,C, C.0,-1 I 

CALL  SYMBOL  (6.32, 1. 17,. 14.5HM/S£C,0.0,5I 

IF  [PLT.EQ.O ) CALL  SYMBOL  (5. 34, .66,. 14, 2HS* ,0.0*2) 

IF  (PLT.EQ.O)  CALL  NUMBER  ( 5 . 62 . . 66, . 14, S , 0. 0 ,- 1 I 
IF  (PLT.EQ.O)  CALL  SYMBOL  ( 6. 32 . . 66, . 14, 5HM/ SEC ,0.0 ,5 1 
IF  (Pl.T.GT.O)  CALL  SYMBOL  ( 5.34, .66,. 14, TH.O.C, 20) 

C •**  PLOT  BRACKETS 

CX  ( 9 ) =»X( N*1 ) 

CX  ( 9 ) -X I N*2 ) 

CY l 9 I *Y  < N* 1 ) 

CY  1 9 ) *Y ( N*2 ) 

00  5 1*1,7 

CX  ( Ii*9X( I l*CX(9)*CX(8) 

5 CY(I)*8Y(I!*CY(9»*CY(8> 

CALL  LINL  (CX.CY, 7,1,0,11 

C ***  PLOT  CURVE 

XMAX*X(N*1)*B.*X(N+2I 
YM AX  Y(NtlI*5.*YlN*2) 

NN*2 

AX(  1>=X(N*1J 
AY  ( l ) -0. J 
AX ( 2 ) *C  ‘ 

AY  (21*0.0 

OR*(Y(N*l)+5.0*Y(N+2))/48. 

IF  ( AX(  2 I *LT  . AX  ( 1 1 ) GO  TO  17 

6 00  7 I *3 c50 

AY( i )»AY( t-11+OR 
IF  (AY( II.GE.YMAX)  GO  TO  8 
AX {!)-(( AY( l ) / A ) **P  »C**P  J **P INV 
IF  (AXUi.Gt.XMAX)  GO  TO  8 
NN*  I 

7 CONTINUE 

8 AX ( NN* 1 I *X ! N* 1 ) 

AX ( NN*2 } *X l N*2 ) 

AY (NN*1 J *YJN*i ) 

AV (NN*2 I “Y( N*2 I 

CALL  LINE  ( AX, AY, NN, 1,0,0) 

C ***  PLOT  ASYMPTOTE 
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NN*  I 

AX(1I*X(N#1> 

AY  1 1 1*Y  ( N*  l ) 

0R»(XMAX-XIN+l))/49. 

00  9 S-2,50 
Axm-Axi  i-n+OR 
IF  ( AXC  n.GE.XHAXJ  GO  TO  10 
AY  1 1 I « AX  I 1 l*A 

IF  IAYII1.GE.YMAX)  GO  TO  10 

9 NN«I 

10  00  11  I-UNN 

AX(II»IAXm-X«N^l)I/XIN*21-.C7 
ATI  I MI  AYm-VlN4l)»/YIM*2)-.035 

11  CALL  SYMBOL  ( AX ( 1 1 , AYIII , . 1 A, IH. , 0.0, 1) 

♦**  PLOT  BOROERS 
CXI6MXIN+1) 

CX|7MX<N*21 
CY  ( 6 ) = YXN«-'l  ) 

CYI7  )=Y(N*21 
OO  12  I*U5 

CXI  I MHDXIIMCXI7MCX16! 

12  CYI  I M0OYI IMCYI71+CYI61 
CALL  LINE  (CX.CY, 5,1,0,11 
CALL  PLOT  1 0. , Y INC ,-3) 

WRITE  (6,20) 

RETURN 

13  SAVE1=X1 11 
SAVE2*Y  <11 
JJ  = 1 

xi  imc 

Y( 1 ) *0« 

GO  TO  1 

14  XI1MSAVE1 
Yl  IMSAVE2 
GO  TO  2 

ASSURES  SPACE  ON  PLOT  BETWEEN  XMIN  AND  X*C 

15  XlN  + lMAA 
YIN»1M0. 

N«N*1 

IF  (XIN).LE.O.)  XINMO. 

CALL  SCALE  IX,8.0,N,ll 
CALL  SCALE  (Y,5.0,N, 1) 

XtNMXIN+1) 

YINMYIN*!! 

XX  N*  l MX  I N+2 1 
YX  N*  1 MY  ‘ N+2  ! 

N»N-l 
GO  TO  3 

16  XIN*i)*PLFXl 
XI N+2 ) *PLOX l 
Y(  N*1  M9, 

YI  N + 2 MPl OY 1 
GO  TO  3 

17  AYI2!*A*IAXI l)**P-C**P)**PINV 
AX  12  MAX  1 1) 

GO  TO  6 

10  K*0 

00  19  J*1*N 

IF  IX! JI.LT.PLFXll  GO  TO  19 
K*K*  1 


non 


AXIKI-XI J1 
AY(XJ-Y(J1 

19  CONTINUE 

axik*i)-x(n«i) 

AX(K+2)-X(N»2) 

AV(K»l)«Y<NM> 

AV(IU2)»Y<N*2> 

CALL  LINE  <AX,AY,K,1,-X,1) 
CO  TO  A 

20  FORMAT  (12H  PLOT  1 OONEI 
END 


SUBROUTINE  PL0T2 
COMMON  (USE  MAIN) 

DIMENSION  DL { 3 i 

DATA  (OL(U. 1-1.3!  /10H , ICK . 10H / 

C ***  PLOT  TITLES. AXES*  ANO  PARAMETERS 

CALL  AXIS  <0.0,0. C,3HC/Xr-3TA.O, 0.0, 0.0, .25) 

CALL  AXIS  (0.0, 0.0, AMY /AX .4, A. 0,90. 0,0.0, .25) 

CALL  SYMBOL  (-.06 , 3. 93 , . 1 A, DL, 0.0 , 30 ) 

CALL  SYMBOL  (4. 07, -. 06,. 14, DL, 90.0, 30) 

CALL  SYMBOL  (-.5,4.5,. 14. TITLE, 0.0,70! 

CALL  SYMBOL  (6.355, 3. 605,. 07, 7HP  1/P, 0.0, 7) 

CALL  SYMBOL  ( 4. 5, 3.5, . 14, TO , 0.0, 22 ) 

CALL  SYMBOL  ( 4. 5 , 2.7, . 14, TA , 0.0, 19) 

CALL  SYMBOL  (4.5, 2.4,.  14,70,0.0,19) 

CALL  SYMBOL  (4.5, 1.9, . 14, 2HP-.0.0, 2) 

CALL  NUMBER  ( 4. 78 , 1. 9, . 1 4,P ,0.0, 1 ) 

IF  IC.EQ, U0R.L.EQ.3I  CALL  SYMBOL  ( 5. 48, 1. 9, • 14  ,TE,0. 0, 71 
CALL  SYMBOL  ( •».  S,  1 .6 , . 14 , 2HA* , 0.0,  2) 

CALL  NUMBER  ( 4 .78, l . 6, . 14, A, 0.0, 2 ) 

IF  (L.E0.2.0R.L.E0.3)  CALL  SYMBOL  1 5. 48, 1. 6,. 14 ,TE ,0.0,71 
CALL  SYMBOL  ( 4 .5,  l. 3 , . 14, 2HC-.0.0, 2 ) 

CALL  NUMBER  (4.78,1.3,.14,C,C.0,-n 
CALL  SYMBOL  ( 5.48 , 1. 3, . 14, 5HM/S6C.O.O, 5 ) 

IF  (PLT.GT.O)  CALL  SYMBOL  ( 4. 50,. 79, . 14, TH ,0.0, 20 1 
C ***  PLOT  CURVE 

AX ( 1 ) -la 
PINV-l./P 
00  l 1-2,50 

1 AX(I)»(1.-FX(I)**P|**PINV 
AX (51) -0.0 

AX (52) >0.0 
AX(53)«. 25 

CALL  LINE  (FX, AX, 51, 1,0,0) 

C ***  PLOT  DATA 

4-.0 

00  2 I-l.N 

IF  (xm.LT.C)  CO  TO  2 
J-J+l 

CXtJ)-C/XtI ) 

CY(J)-Y(  l)/CA*X(I)) 

2 CONTINUE 
CXU  + ll-0.0 
CX  (J  + 2I-.25 
CYU  + U-0.0 
CVI4»2)-.25 


48 


nn  o 


CALL  LINE  ICX.CY, J, 1,-1, l) 
c PLOT  80R0ERS 

CX(6|«0.0 
CX (7) *.25 
CY(6) *0.0 
CY 1 7 1 *»  25 
DO  3 1-1,5 

CXCI)=BOXC ! I *CX ( 7 I *CXI M 

3 CYm*B0im*CY|7|*CY(6l 
CALL  LINE  ICX.CY, 5.1» 0,l» 
CALL  PLOT  IO..YINC.-3) 
WRITE  (6,4) 

RETURN 

C 

4 FORMAT  ( 12H  PLOT  2 DONE) 
END 


SUBROUTINE  PLOT3 
COMMON  (USE  MAINI 
OIKENSION  DYl 100) 

DO  1 1*1, N 

1 OY  l 1 > —X  ( ll-Yd) 

IF  IF  tXAX.GT ,0)  GO  TO  10 
XIN*2I=C 
OY ( N + 2 ) =C 
X l N* 1 ( =0 . 

OY  ( N » 1 J * 3 » 

N-N+2 

CALL  SCALE  (X.6.0.N, l) 

CALL  SCALE  { OY , 4*0 ,N . l I 
XCN-l)=X(N*l) 

XI N J *XI N+2 I 
OYIN-lI-OYlN+l) 

DY I N) =DY ( N*2 ) 

N*N-2 

C *•*  PLOT  TITLES, AXES,  AND  P.  URETERS 

2 CALL  AXiS  1 0 .0 ,0 .0 , 9HX  I M/SEC 5 ,-9, 6.0 ,0. 0, X I N* l ) , X ( N+2 ) ) 
CALC  AXIS  l0.0,0.0,TG,ll,4.j,90.0,0Y<N*l>,DY(N*2)l 

CALL  SYMBOL  (-.5, 5.5,. 14, TITLE, C. 0,701 
CALL  SYMBOL  ( 1.98,4.72,. 14, 4HX-Y=, 0.0, 4) 

CALL  SYMBOL  ! 2 . 82 , 4. 79 , . 1 4 , OHX , 0 X C, 0.0, 81 
CALL  SYMBOL  (3,41,4.86,.07,5H<  <,0.0, 5) 

CALL  SYMBOL  ( 3. 4 1 , 4. 79 , . 07 , 5H  = =,C.0,5I 

CALL  SYMBOL  ( 3 .59, 4. 685, .07, TO, C. 0, 20  I 
CALL  SYMBOL  ( 2 . 82 , 4, 58 , . 14, TF , 0. 0, 20 ) 

CALL  SYMBOL  (5,52,3. 0,. 14, TA.C.O, 19) 

CALL  SYMBOL  ( 5 . 52 , 2. 8 , . 1 4, TC, 0. V, 1 9 1 
CALL  SYMBOL  ( 6 . 5 • 2 . 3 , . 1 4, 2HP  = , 0 .0, 2 ) 

CALL  NUMBER  (6. 78, 2. 3, . 14, P, 0.0, II 

IF  I L • EQ • 1 • OR • L . EQ • 3 i CALL  SYMBOL  i 7. 48 , 2* 3 , » 14 , It , 0* 0, I ) 
CALL  SYMBOL  ( 6. 5 , 2.0 , . 14 , 2HA  = ,C,0» 2 1 
CALL  NUMBER  ( 6.78, 2.0, . 14, A, 0.0, 21 

IF  IL.E0.2.JR.L.E0.3)  CALL  SYMBOL  17. 48 ,2. 0, .14 ,TE ,0,0,71 
CALL  SYMBOL  ( 6.5, 1.7,, 14, 2HC*. 0,0, 21 
CALL  NUMBER  (6.78,  1.7, ,14,C,0.G,-l» 

CALL  SYMBOL  (7.34, 1.7,.14,6HM/SEC  ,0.0,61 
IF  IPLT.GT.O)  CALL  SYMBOL  I 6. 5, 1. 19, . 14, TH.0,0 , 20) 

C ***  PLOT  BRACKETS 
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CX(8)-X(N«l) 

CXm*X(N»2) 

cvm-ov<N+n 

CYI9)-0Y(N*2> 

00  3 I>1,7 

CXni-'BX(I)-3.08)*CXOt*CX(6l 

3 CY(l)«(BY<  I )«2.54i*C.Ym«CYtBt 
CALL  LINE  ( CX,CY, 7, 1,0, l ) 

C • ••  PLOT  CURVE 

NN-2 
CXID-O. 

cvm-o.  | 

cxm-c  . ; 

CVI2I-C  I 

RHAX»4.*CV(N*2) 

2HAX*6»*X(N+2)  i 

02*UNAX-CX(2))/43.  . j 

PINV-l./P 
10  4 1-3,50 

cxi  n-cxi  mud: 

IF  (CX(I).GT.2HAX}  GO  TO  5 

cyu)-cxi  i»-a*icx(  ij**p-c**pi**pinv  j 

IF  ( CY ( I } .GT.RMAX  9 GO  TO  5 

NN-1  ! 

4 CONTINUE  1 

5 CXINN+D-X1N+H  j 

CX(NN+2)*X(N+2) 

CY(NN*l)*OY(N*l)  1 

CY  (NN+2I-0YIN+2) 

CALL  LINE  (CX.CY.NN, 1,0,01 
C PLOT  ASYMPTOTE 

0-(l.-A|*2HAX 
IF  (O.GT.RMAXJ  GO  TO  11 

6 2- 2 MAX/ FLO AT (NN-1 ) 

cxm-o. 

CYUI-0. 

DO  7 1*2, NN 
CXUi-CXI  I-l)*2 

7 CYU  >-(  l.-A)*CX(I) 

00  8 1*1 .NN 

CX(I)-(CX( I )-X(N+l))/XSN+2»-.07 
CY ( 1 1* (CY( I »-DY(N*l) }/0Y(N*2l-.035 
CALL  SYMBOL  <CX( I ) ,CY( 1 1 , . U, 1H, ,0.0, 11 

8 CONTINUE 

C •**  PLOT  DATA 

CALL  LINE  ( X , OY ,N, 1,-1, II 
C ***  PLOT  BORDERS 

CX(6t*X(N«l) 

CX ( 7 > *x; N*2  J 
CV(6)»Y(NUi 
CY(7)«Y(N*2) 

DO  9 1*1,5 

CXm»BCXCn«CX(7)+CX(6) 

9 CY ( 1 1* BOY  1 1 } *CY ( 7 I +CY( 6 1 
CALL  LINE  (CX.CY, 5,1, 0,11 
YORIG— 2,*YINC 

CALL  PLOT  (XINC.YORIG.-35 
WRITE  (6,12) 

WRITE  (6,13) 

RETURN 
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RETURN 

10  XtN+ll-0. 

XIN+2J-PL0X3 

0VtN*l>»0. 

OV  CN* 2) ■PLOY 3 
GO  TO  2 

11  ZMAX»RMAX/( l.-At 
GO  TO  6 

12  FORMAT  ( 12H  PLOT  3 DONE} 

13  FORMAT  < 1HU 
END 


